# Analysis commands for:
# ``The Conditional Impact of Military Intervention on Internal Armed Conflict Outcomes
# Patricia L. Sullivan and Johannes Karreth
# Conflict Management and Peace Science 32(3): 269-288.

# Last updated: November 11, 2015 (added info to Table 3)
# Data and analysis commands stored at:
# https://dataverse.harvard.edu/dataverse/jkarreth

##############################
# Setup
##############################

# setwd("...")
# sink(file = "cimi.log.txt", append = FALSE, type = "output") # For log file
dat <- read.csv("cimi.data.csv")

# Note: Coding of outcome_d (d = (conflict-)dyad)
# 0: low activity
# 1: rebel victory
# 2: settlement
# 3: govt victory

##############################
# Table 1
##############################

library(nnet)

dat.m1 <- dat[is.na(dat$outcome_d) == FALSE & is.na(dat$intervention) == FALSE, c("outcome_d", "intervention")]
m1 <- multinom(outcome_d ~ intervention, data = dat.m1)

dat.m2 <- dat[is.na(dat$outcome_d) == FALSE & is.na(dat$intervention) == FALSE & is.na(dat$fightcaphigh) == FALSE & is.na(dat$mi_fightcap) == FALSE, ]
m2 <- multinom(outcome_d ~ intervention + fightcaphigh + mi_fightcap, data = dat.m2)

dat.m3 <- dat[is.na(dat$outcome_d) == FALSE & is.na(dat$intervention) == FALSE & is.na(dat$fightcaphigh) == FALSE & is.na(dat$mi_fightcap) == FALSE & is.na(dat$troops2rebs) == FALSE & is.na(dat$coup) == FALSE & is.na(dat$land) == FALSE & is.na(dat$lngdp) == FALSE & is.na(dat$lnyears) == FALSE & is.na(dat$postCW) == FALSE, ]
m3 <- multinom(outcome_d ~ intervention + fightcaphigh + mi_fightcap + troops2rebs + coup + land + lngdp + lnyears + postCW, data = dat)

library(memisc)
mtable(m1, m2, m3, digits = 3)

# Model checks

# dat.m3.noint <- dat[is.na(dat$outcome_d) == FALSE & is.na(dat$intervention) == FALSE & is.na(dat$fightcaphigh) == FALSE & is.na(dat$troops2rebs) == FALSE & is.na(dat$coup) == FALSE & is.na(dat$land) == FALSE & is.na(log(dat$lngdp)) == FALSE & is.na(dat$lnyears) == FALSE & is.na(dat$postCW) == FALSE, ]
# m3.noint <- multinom(outcome_d ~ intervention + fightcaphigh + troops2rebs + coup + land + lngdp + lnyears + postCW, data = dat.m3.noint)
# 
# library(separationplot)
# sp.categorical(pred = fitted(m1), actual = dat.m1$outcome_d)
# sp.categorical(pred = fitted(m2), actual = dat.m2$outcome_d)
# sp.categorical(pred = fitted(m3), actual = dat.m3$outcome_d)
# sp.categorical(pred = fitted(m3.noint), actual = dat.m3.noint$outcome_d)
# 
# anova(m3, m3.noint)

# PRE for each model by hand:

# mode <- function(x) {
#   ux <- unique(x)
#   ux[which.max(tabulate(match(x, ux)))]
# }

predicted.m1 <- predict(m1)
correct.m1 <- ifelse(predicted.m1 == dat.m1$outcome_d, 1, 0)
table(correct.m1)
modal.m1 <- max(table(dat.m1$outcome_d))
corr.pred.m1 <- sum(correct.m1) # / length(dat.m1$outcome_d)
pre.m1 <- 100 * ((corr.pred.m1 - modal.m1) / (nrow(dat.m1) - modal.m1))
print(pre.m1, digits = 3)

predicted.m2 <- predict(m2)
correct.m2 <- ifelse(predicted.m2 == dat.m2$outcome_d, 1, 0)
table(correct.m2)
modal.m2 <- max(table(dat.m2$outcome_d))
corr.pred.m2 <- sum(correct.m2) # / length(dat.m2$outcome_d)
pre.m2 <- 100 * ((corr.pred.m2 - modal.m2) / (nrow(dat.m2) - modal.m2))
print(pre.m2, digits = 3)

predicted.m3 <- predict(m3)
correct.m3 <- ifelse(predicted.m3 == dat.m3$outcome_d, 1, 0)
table(correct.m3)
modal.m3 <- max(table(dat.m3$outcome_d))
corr.pred.m3 <- sum(correct.m3) # / length(dat.m3$outcome_d)
pre.m3 <- 100 * ((corr.pred.m3 - modal.m3) / (nrow(dat.m3) - modal.m3))
print(pre.m3, digits = 3)

# Effects plot (Model 3)

# dat.m3.f <- dat.m3
# dat.m3.f$GovtIntervention <- as.factor(dat.m3$intervention)
# levels(dat.m3.f$GovtIntervention) <- c("No", "Yes")
# dat.m3.f$RebelStrength <- as.factor(dat.m3$fightcaphigh)
# levels(dat.m3.f$RebelStrength) <- c("Weak", "Strong")
# dat.m3.f$RebelIntervention <- as.factor(dat.m3$troops2rebs)
# levels(dat.m3.f$RebelIntervention) <- c("No", "Yes")
# dat.m3.f$Coup <- as.factor(dat.m3$coup)
# levels(dat.m3.f$Coup) <- c("No", "Yes")
# dat.m3.f$Separatist <- as.factor(dat.m3$land)
# levels(dat.m3.f$Separatist) <- c("No", "Yes")
# dat.m3.f$PostCW <- as.factor(dat.m3$postCW)
# levels(dat.m3.f$PostCW) <- c("No", "Yes")
# dat.m3.f$Outcome <- as.factor(dat.m3$outcome_d)
# levels(dat.m3.f$Outcome) <- c("Low activity", "Rebel victory", "Negotiated settlement", "Gov't victory")
# 
# m3.f <- multinom(Outcome ~ GovtIntervention*RebelStrength + RebelIntervention + Coup + Separatist + lngdp + lnyears + PostCW, data = dat.m3.f)

# library(effects)
# 
# plot(effect("Intervention", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("RebelStrength", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("GovtIntervention*RebelStrength", m3.f), typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("RebelIntervention", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("Coup", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("Separatist", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("lngdp", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("lnyears", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)
# plot(effect("PostCW", m3.f), levels = 90, typical = "median", levels = 90, as.table = TRUE, rescale.axis = FALSE)


##############################
# Table 2
##############################

# Zelig setup

library(Zelig)
library(ZeligChoice)
dat$outcomeZfoo <- as.numeric(dat$outcome_d)
dat$outcome_d.Z <- ifelse(dat$outcomeZfoo == 0, 4, dat$outcomeZfoo)
dat$outzomeZfoo <- NULL

dat.m1.Z <- dat[is.na(dat$outcome_d.Z) == FALSE & is.na(dat$intervention) == FALSE, c("outcome_d.Z", "intervention")]
dat.m2.Z <- dat[is.na(dat$outcome_d.Z) == FALSE & is.na(dat$intervention) == FALSE & is.na(dat$fightcaphigh) == FALSE & is.na(dat$mi_fightcap) == FALSE, c("outcome_d.Z", "intervention", "fightcaphigh", "mi_fightcap")]
dat.m3.Z <- dat[is.na(dat$outcome_d.Z) == FALSE & is.na(dat$intervention) == FALSE & is.na(dat$fightcaphigh) == FALSE & is.na(dat$mi_fightcap) == FALSE & is.na(dat$troops2rebs) == FALSE & is.na(dat$coup) == FALSE & is.na(dat$land) == FALSE & is.na(dat$lngdp) == FALSE & is.na(dat$lnyears) == FALSE & is.na(dat$postCW) == FALSE, c("outcome_d.Z", "intervention", "fightcaphigh", "mi_fightcap", "troops2rebs", "coup", "land", "lngdp", "lnyears", "postCW")]

# Note: Coding of outcome_d.Z
# 1: rebel victory
# 2: settlement
# 3: govt victory
# 4: low activity

m1.out <- zelig(as.factor(outcome_d.Z) ~ intervention, data = dat.m1.Z, model = "mlogit")
m2.out <- zelig(as.factor(outcome_d.Z) ~ intervention + fightcaphigh + mi_fightcap, data = dat.m2.Z, model = "mlogit")
m3.out <- zelig(as.factor(outcome_d.Z) ~ intervention + fightcaphigh + mi_fightcap + troops2rebs + coup + land + lngdp + lnyears + postCW, data = dat.m3.Z, model = "mlogit")

# Baseline probabilities

baseline.x <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
baseline.sim <- sim(m3.out, x = baseline.x)
# summary(baseline.sim)

# First differences

# Pro-govt intervention when rebels are weak

miweakrebs.lo <- setx(m3.out, intervention = 0, fightcaphigh = 0, mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
miweakrebs.hi <- setx(m3.out, intervention = 1, fightcaphigh = 0, mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
miweakrebs.sim <- sim(m3.out, x = miweakrebs.lo, x1 = miweakrebs.hi)      
# summary(miweakrebs.sim)

miweakrebs.fd <- apply(miweakrebs.sim$qi$fd, 2, mean)
miweakrebs.lower <- apply(miweakrebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
miweakrebs.upper <- apply(miweakrebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))

# Strong rebels, no pro-govt intervention

strongrebs.lo <- setx(m3.out, intervention = 0, fightcaphigh = 0, mi_fightcap = 0, troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
strongrebs.hi <- setx(m3.out, intervention = 0, fightcaphigh = 1, mi_fightcap = 0, troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
strongrebs.sim <- sim(m3.out, x = strongrebs.lo, x1 = strongrebs.hi) 
# summary(strongrebs.sim)
strongrebs.fd <- apply(strongrebs.sim$qi$fd, 2, mean)
strongrebs.lower <- apply(strongrebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
strongrebs.upper <- apply(strongrebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))

# Pro-govt intervention when rebels are strong

mistrongrebs.lo <- setx(m3.out, intervention = 0, fightcaphigh = 1, mi_fightcap = 0, troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
mistrongrebs.hi <- setx(m3.out, intervention = 1, fightcaphigh = 1, mi_fightcap = 1, troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
mistrongrebs.sim <- sim(m3.out, x = mistrongrebs.lo, x1 = mistrongrebs.hi) 
# summary(mistrongrebs.sim)
mistrongrebs.fd <- apply(mistrongrebs.sim$qi$fd, 2, mean)
mistrongrebs.lower <- apply(mistrongrebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
mistrongrebs.upper <- apply(mistrongrebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))

# Pro-rebel intervention

troops2rebs.lo <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = 0, coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
troops2rebs.hi <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = 1, coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
troops2rebs.sim <- sim(m3.out, x = troops2rebs.lo, x1 = troops2rebs.hi)      
# summary(troops2rebs.sim)
troops2rebs.fd <- apply(troops2rebs.sim$qi$fd, 2, mean)
troops2rebs.lower <- apply(troops2rebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
troops2rebs.upper <- apply(troops2rebs.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))

# Coup

coup.lo <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = 0, land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
coup.hi <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = 1, land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
coup.sim <- sim(m3.out, x = coup.lo, x1 = coup.hi)      
# summary(coup.sim)
coup.fd <- apply(coup.sim$qi$fd, 2, mean)
coup.lower <- apply(coup.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
coup.upper <- apply(coup.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))

# Land

land.lo <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = 0, lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
land.hi <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = 1, lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = median(dat.m3.Z$postCW))
land.sim <- sim(m3.out, x = land.lo, x1 = land.hi)      
# summary(land.sim)
land.fd <- apply(land.sim$qi$fd, 2, mean)
land.lower <- apply(land.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
land.upper <- apply(land.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))


# Duration

dur.lo <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = quantile(dat.m3.Z$lnyears, probs = c(0.25)), postCW = median(dat.m3.Z$postCW))
dur.hi <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = quantile(dat.m3.Z$lnyears, probs = c(0.75)), postCW = median(dat.m3.Z$postCW))
dur.sim <- sim(m3.out, x = dur.lo, x1 = dur.hi)
# summary(dur.sim)
dur.fd <- apply(dur.sim$qi$fd, 2, mean)
dur.lower <- apply(dur.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
dur.upper <- apply(dur.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))


# Post-Cold War

postcw.lo <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = 0)
postcw.hi <- setx(m3.out, intervention = median(dat.m3.Z$intervention), fightcaphigh = median(dat.m3.Z$fightcaphigh), mi_fightcap = median(dat.m3.Z$mi_fightcap), troops2rebs = median(dat.m3.Z$troops2rebs), coup = median(dat.m3.Z$coup), land = median(dat.m3.Z$land), lngdp = median(dat.m3.Z$lngdp), lnyears = median(dat.m3.Z$lnyears), postCW = 1)
postcw.sim <- sim(m3.out, x = postcw.lo, x1 = postcw.hi)
# summary(postcw.sim)
postcw.fd <- apply(postcw.sim$qi$fd, 2, mean)
postcw.lower <- apply(postcw.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.05)))
postcw.upper <- apply(postcw.sim$qi$fd, 2, function(x) quantile(x, probs = c(0.95)))


fd.dat <- rbind(t(baseline.sim$stats[[1]])[c(1, 4, 5), ], miweakrebs.fd, miweakrebs.lower, miweakrebs.upper, strongrebs.fd, strongrebs.lower, strongrebs.upper, mistrongrebs.fd, mistrongrebs.lower, mistrongrebs.upper, troops2rebs.fd, troops2rebs.lower, troops2rebs.upper, coup.fd, coup.lower, coup.upper, land.fd, land.lower, land.upper, dur.fd, dur.lower, dur.upper, postcw.fd, postcw.lower, postcw.upper)
colnames(fd.dat) <- c("Rebel victory", "Settlement", "Gov't victory", "Low activity")

# Rounded first differences & CIs for Table 2

fd.dat.rounded <- apply(fd.dat, c(1, 2), function(x) round(x, digits = 2))
library(reshape2)
table2 <- melt(fd.dat.rounded)
is.wholenumber <-
  function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol
table2$rownum <- c(1:108)
Lower95 <- table2[is.wholenumber((table2$rownum + 1) / 3) == TRUE, c(1:3)]
Upper95 <- table2[is.wholenumber((table2$rownum) / 3) == TRUE, c(1:3)]
table2 <- table2[is.wholenumber((table2$rownum - 1) / 3) == TRUE, c(1:3)]
table2$Lower95 <- Lower95$value
table2$Upper95 <- Upper95$value
table2$Variable <- rep(c("Baseline", "Pro-Govt Intervention vs Weak Rebels", "Strong Rebels", "Pro-Govt Intervention vs Strong Rebels", "Pro-Rebel Intervention", "Coup", "Secessionist Conflict", "Duration (25th to 75th percentile", "Post-Cold War"), 4)
table2$Outcome <- table2$Var2
table2$Mean <- table2$value
table2 <- table2[, c(6, 7, 8, 4, 5)]
table2

##############################
# Table 3 (unit of analysis: conflict-year)
##############################

# Note: Coding of outcome_dy (dy = (conflict-)dyad-year)
# 0: continue
# 1: rebel victory
# 2: low activity
# 3: settlement
# 4: govt victory


m3.dy <- multinom(outcome_dy ~ exgov_combat_dy + strongrebs_dy + mi_strongrebs_dy + exreb_combat_dy + coup_dy + land_dy + lngdp_dy + lnyears_dy + postCW_dy, data = dat)

mtable(m3.dy)

sessionInfo()
sink()